import cv2
from matplotlib import pyplot as plt
import numpy as np
import os
import scipy.signal as sgn
# Ivo's implemetation of Ellipse class
class Ellipse:
C : np.array #homogeneous conic matrix x'*C*x=0
def __init__( self, C : np.array ):
self.C = C
@classmethod
def estimate_ellipsoid( cls, points : np.array ):
"""
C = estimate_ellipsoid_2D ( points )
----------------------------------
estimates a homogeneous ellipsoid C
x'*C*x = 0
from a number of point samples
(in a least squares sense)
format of points
[ x^1_1 . . . . . x^N_1 ]
[ . . ]
points = [ . . ]
[ . . ]
[ x^1_n . . . . . x^N_n ]
[ 1 . . . . . 1 ]
where each point is element R^n
and there are N points, the points
have to be in homogeneous coordinates.
"""
n = points.shape[0]
print(n)
num_entries = int((n+1) * n/2)
if points.shape[1] < num_entries-1:
print('insufficient number of points, need at least {}'.format( num_entries-1 ) )
# -------------- estimate ellipse ---------------
# assemble linear system
A = np.zeros( [points.shape[ 1 ], num_entries ] )
#for all points
for i in range( 0, points.shape[ 1 ] ):
#extract point
x = np.array([points[:, i]]).T
#print(x)
X = np.matmul( x, x.T )
#print(X)
X2 = X.T + X - np.diag( np.diag( X ) )
#print(X2)
row = X2[ np.where( np.triu( np.ones( X2.shape ) ) ) ]
#print(row)
A[ i, : ] = row
#print(A)
#solve for null space
U,D,V = np.linalg.svd( A )
ell = V[-1,:]
#create conic matrix
Cd = np.zeros( X2.shape )
Cd[ np.where( np.triu( np.ones( X2.shape ) ) ) ] = ell
C = Cd + Cd.T - np.diag( np.diag ( Cd ) )
#detC = np.linalg.det( C[0:-1,0:-1] )
#print(detC)
if Ellipse.is_hyperboloid( C ):
print('Warning: result is a hyperboloid!')
if Ellipse.is_paraboloid( C ):
print('Warning: result is a paraboloid!')
return C
@classmethod
def cov_mean_to_ellipsoid( cls, cov, mean, k ):
"""
C = cov_mean_to_ellipsoid_homog ( Cov, Mean, k )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from covariance matrix and mean value
of a gaussian distribution
k is a confidence factor determibed by the chi^2
distribution. The probability that x is inside
the ellipsoid is P_{chi^2}(k,2). Some values
for k are 2.41 (70%) and k=5.99 (95%).
cov : np.array - covariance matrix (non-homogeneous)
mean : np.array - mean value (vector)
k : float - confidence factor, see above
"""
#column vector, assume we get a 1D array
mean = np.array( [mean] ).T
SIGMAinv = np.linalg.inv( cov )
lastentry = np.matmul( mean.T, np.matmul( SIGMAinv, mean ) ) - k
C = np.vstack( [ np.hstack( [ SIGMAinv, -np.matmul(SIGMAinv, mean)] ),
np.hstack( [ -np.matmul(mean.T,SIGMAinv), lastentry] ) ] )
return C
@classmethod
def srt_to_ellipsoid( cls, S, R, T ):
"""
C = srt_to_ellipsoid_homog ( S, R, T )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from translation, rotation and scale parts as
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
see: Ivo Ihrke: "Some Notes on Ellipses", TechReport, 2004
"""
#function C = srt_to_ellipsoid_homog ( S, R, T )
invT = np.linalg.inv( T )
invR = R.T
invS = np.linalg.inv( S )
left = np.matmul( invT.T, np.matmul( invR.T, invS.T ) )
right = np.matmul( invS, np.matmul( invR, invT ) )
C = np.matmul( left, right )
return C
@classmethod
def decompose_ellipsoid( cls, C ):
"""
[S,R,T] = decompose_ellipsoid_homog ( C )
----------------------------------
decomposes a homogeneous ellipsoid C
x'*C*x = 0
into translation, rotation and scale parts
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
Note that C \neq (T^-1 * R^-1 * S^-1) ( S * R * T )
in some cases because the orientation of
the ellipsoid's major axis is not fixed.
The function generates a right handed system.
"""
#print(C)
C_ih = C[0:-1,0:-1]
#print(C_ih)
t = -np.matmul( np.linalg.pinv( C_ih ), C[0:-1,-1] )
T = np.eye( C.shape[ 0 ] )
T[0:-1, -1 ] = t
U,D,V = np.linalg.svd( C_ih )
R = np.eye( C.shape[ 0 ] )
R[0:-1,0:-1] = V
#S = R' * T' * C * T * R;
S = np.matmul( R.T, np.matmul( T.T, np.matmul( C, np.matmul( T, R ) ) ) )
S = S / S[-1,-1]
S = np.sqrt( np.diag(1/abs(np.diag(S))))
return S, R, T
@classmethod
def is_ellipsoid( cls, C ):
"""
bool = is_ellipsoid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is an ellipsoid
"""
#function val = conic_is_ellipsoid_homog ( C )
detC = np.linalg.det( C[0:-1,0:-1] )
return detC > 0
@classmethod
def is_hyperboloid( cls, C ):
"""
bool = is_hyperboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a hyperboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
return detC < 0
@classmethod
def is_paraboloid( cls, C ):
"""
bool = is_paraboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a paraboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
#maybe with EPS ?
return detC == 0
def get_center( self ):
S,R,T = self.decompose_ellipsoid( self.C )
return [ T[2,0], T[2,1] ]
def plot( self, ax, **kwargs ):
"""
plot_ellipsoid_2D ( C, style, N )
----------------------------------
plots a homogeneous ellipsoid C
x'*C*x = 0
style is a matlab plot style, e.g. 'b-' (default)
N is the number of polygon segments to
approximate the ellipse
"""
#function plot_ellipsoid_homog_2d ( C, style, N )
if not 'style' in locals():
style = 'b-'
if not 'N' in locals():
N = 100
angles = np.linspace( -np.pi, np.pi, N )
# generate a unit circle or a unit hyperbola
if ( Ellipse.is_ellipsoid( self.C ) ):
#uc = [ cos(x); sin(x); ones(size(x)) ];
x = np.cos( angles )
y = np.sin( angles )
w = np.ones( angles.shape )
uc = np.vstack([x,y,w])
else:
EPS = 0.1
x = np.linspace( -np.pi/2+EPS, np.pi/2-EPS, int( N/2 ) )#-(pi/2-EPS):((pi-2*EPS)/(N/2)):(pi/2-EPS);
x2 = np.linspace( np.pi/2+EPS, 3*np.pi/2-EPS, int( N/2 ) ) #x2 = (pi/2+EPS):((pi-2*EPS)/(N/2)):(3*pi/2-EPS);
uc = np.vstack( [ 1/np.cos(x), np.tan(x), np.ones( x.shape ) ] )
uc2 = np.vstack( [ 1/np.cos(x2), np.tan(x2), np.ones( x2.shape ) ] )
#decompose
S,R,T = Ellipse.decompose_ellipsoid( self.C )
#transform unit circle or parabola
transform = np.matmul( T, np.matmul( R, S ) )
ell = np.matmul( transform, uc )
if 'uc2' in locals():
ell2 = np.matmul( transform, uc2 )
ax.plot( ell2[0,:], ell2[1,:], style )
#plot conic
ax.plot( ell[0,:], ell[1,:], style )
#plot coordinate system
tx = T[0,2]
ty = T[1,2]
# print(tx, ty)
smin = S[0,0]
smaj = S[1,1]
ax.plot( np.array( [ tx, tx+smin*R[0,0] ] ), np.array( [ ty, ty+smin*R[1,0] ] ), 'r' )
ax.plot( np.array( [ tx, tx+smaj*R[0,1] ] ), np.array( [ ty, ty+smaj*R[1,1] ] ), 'b' )
ax.set_aspect('equal', adjustable='box')
return (tx, ty)
def fit_ellipse(path_dir_w_time, ifHull=False, date_now=""):
""" Find and fit an ellipse in an image.
it should be run on a normalized image I (normalization should exclude dead pixels which are artificially maximum)
Parameters
----------
path_dir_w_time : str
_description_
date_now : str
_description_
Returns
-------
Ellipse
Ellipse in an image.
"""
I_temp = plt.imread(path_dir_w_time)[:,:,1].astype(float)/255
I = old_rplc_dead_px_med(I_temp)
ellipsoid_points=[]
if ifHull == False:
#every Nth scanline should be sampled
START_SCANLINE = 50 #there are some weird effects on the zero scanline that lead to detection of erroneous points
DELTA_SCANLINE = 100
#need to filter enough such that the big jump into and out of the disc has indeed the largest magnitude derivative
FILTERWIDTH=250
x=np.linspace(-4,4,FILTERWIDTH)
gauss=np.exp(-x**2)
plt.figure(figsize=[18,12])
for sl in range( START_SCANLINE, I.shape[0]-START_SCANLINE, DELTA_SCANLINE ):
#print(f'Checking Scanline {sl} ...' )
#line = self.I[sl]
line = I[sl, :]
filtered_line = sgn.convolve(line,gauss,mode='same')
#crude filter for line intersecting the ring structure
if np.amax(filtered_line) > 0.9:
#print("- scanline intersection ")
filtered_deriv = np.diff( filtered_line )
#crude extraction - this will probably need improvement
max_pos = np.argmax( filtered_deriv )
min_pos = np.argmin( filtered_deriv )
ellipsoid_points.append([max_pos,sl])
ellipsoid_points.append([min_pos,sl]) #add two points on perimeter
else:
#print("- no intersection ")
pass
print(ellipsoid_points)
else:
img = np.array(I*255, dtype=np.uint8)
cnt_img, final_cnt = contour_points(img)
if final_cnt is None:
return None, None
elif final_cnt.shape[0] < 5:
return None, None
else:
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
# print(ellipsoid_points)
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
print(ellipsoid_points)
#plt.clf()
#plt.imshow(self.I)
#plt.show()
for i in range( 0, ellipsoid_points.shape[1] ):
plt.plot(ellipsoid_points[0,i],ellipsoid_points[1,i],'rx')
E = Ellipse( Ellipse.estimate_ellipsoid( ellipsoid_points ) )
S,R,T=Ellipse.decompose_ellipsoid(E.C)
plt.plot(T[0,2],T[1,2],'gx')
plt.imshow(I)
center = E.plot(plt.gca())
plt.title(f'ellipsoide fit {1}')
# plt.savefig(f"{path_dir_w_time}/{date_now}_ledpos_{self.x_pos}_{self.y_pos}_{self.channel}_ellipsoide_fit.png")
return E, center
def old_rplc_dead_px_med(I_temp):
"""Replace dead pixels by median
Parameters
----------
I_temp : ndarray
image to be processed
Returns
-------
ndarray
clean image where dead pixels are replaced by median
"""
Imed = sgn.medfilt2d( I_temp, kernel_size=5 )
maxI = np.amax( Imed ) #normalize robustly
I_temp = I_temp / maxI
#replace dead pixels by median filtered values
I_temp[ I_temp > 1 ] = Imed[ I_temp > 1 ]
return I_temp
def contour_points(g_img):
kernel = np.ones((3, 3), np.uint8)
gradient = cv2.morphologyEx(g_img, cv2.MORPH_GRADIENT, kernel)
blur = cv2.GaussianBlur(gradient, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
gradient_img = np.zeros_like(g_img, dtype=np.uint8)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 2:
return None, None
cv2.drawContours(gradient_img, contours, -1, 255, 2)
kernel = np.ones((5, 5), np.uint8)
cv2.dilate(gradient_img, kernel, iterations = 1)
hull_img = np.zeros_like(g_img, dtype=np.uint8)
blur = cv2.GaussianBlur(g_img, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 2:
return None, None
area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
max_cntr_index = np.argmax(area_cntrs)
cnt = contours[max_cntr_index]
hull = cv2.convexHull(cnt)
cv2.drawContours(hull_img, [hull], 0, 125, 2)
op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
op[:10, :] = 0
op[-10:, :] = 0
kernel = np.ones((3, 3), np.uint8)
opening = cv2.morphologyEx(op, cv2.MORPH_OPEN, kernel)
x, y = np.where(opening > 1)
final_cnt = np.vstack([y, x]).T
# final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=int(final_cnt.shape[0]/10)))), axis=0)
imgs = np.hstack([gradient_img, hull_img, opening])
if final_cnt.shape[0] > 4:
return opening, final_cnt
else:
return None, None
dir_path = "Z:/CSE/CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_08_25/2023_08_25_09_31_25/"
file_path = "2023_08_25_09_31_25_img_iso_100_shutter_2_x_16_y_16_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
imgs, cnt = contour_points(g_img)
kernel = np.ones((25, 25),np.uint8)
cv2.dilate(imgs, kernel, iterations = 1)
fig, ax = plt.subplots(1, 1, figsize=(15, 15))
plt.imshow(imgs, cmap='gray')
plt.show()
ellipse = cv2.fitEllipse(cnt)
cv2.ellipse(g_img, ellipse, 255, 2)
fig, ax = plt.subplots(1, 1, figsize=(15, 15))
ax.imshow(g_img, cmap='gray')
plt.show()
print(ellipse)
((2971.643798828125, 913.9298095703125), (5347.255859375, 5739.8564453125), 169.41822814941406)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff"
# dir_path = "Z:/CSE/CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_08_25/2023_08_25_09_31_25/"
# file_path = "2023_08_25_09_31_25_img_iso_100_shutter_2_x_16_y_16_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
g_img
array([[ 0, 0, 0, ..., 0, 0, 0],
[ 0, 0, 0, ..., 0, 0, 0],
[ 0, 0, 0, ..., 0, 0, 0],
...,
[ 887, 153, 1258, ..., 0, 0, 0],
[1531, 1343, 2342, ..., 0, 0, 0],
[1369, 1267, 1454, ..., 0, 0, 0]], dtype=uint16)
plt.imshow(g_img, cmap='gray')
<matplotlib.image.AxesImage at 0x208ad516f10>
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
plt.imshow(g_img_8, cmap='gray')
<matplotlib.image.AxesImage at 0x208b074b070>
kernel = np.ones((15, 15), np.uint8)
dilation = cv2.dilate(g_img_8, kernel, iterations = 1)
plt.imshow(dilation, cmap='gray')
plt.show()
kernel = np.ones((15, 15), np.uint8)
erosion = cv2.erode(g_img_8, kernel, iterations = 1)
plt.imshow(erosion, cmap='gray')
plt.show()
plt.imshow((np.fabs(-erosion.astype(np.int16) + dilation.astype(np.int16))).astype(np.uint8), cmap='gray')
plt.show()
kernel = np.ones((15, 15), np.uint8)
gradient = cv2.morphologyEx(g_img_8, cv2.MORPH_GRADIENT, kernel)
# g_img_2 = g_img.cop
plt.imshow(gradient, cmap='gray')
<matplotlib.image.AxesImage at 0x208b23ed430>
# kernel = np.ones((25, 25), np.uint8)
# opening = cv2.morphologyEx(gradient, cv2.MORPH_OPEN, kernel)
# plt.imshow(opening, cmap='gray')
blur = cv2.GaussianBlur(gradient+20, (5, 5), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
gradient_img = np.zeros_like(g_img, dtype=np.uint8)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 2:
print("Errors")
cv2.drawContours(gradient_img, contours, -1, 255, 2)
kernel = np.ones((35, 35),np.uint8)
gradient_img = cv2.dilate(gradient_img, kernel, iterations = 1)
plt.imshow(gradient_img, cmap='gray')
<matplotlib.image.AxesImage at 0x208d1f2adf0>
ch_img_1 = np.zeros_like(img, dtype=np.uint8)
ch_img_1[:, :, 0] = gradient_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_8+50
new_img = np.array((ch_img_1 + ch_img_2), dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
hull_img = np.zeros_like(g_img_8, dtype=np.uint8)
blur = cv2.GaussianBlur(g_img_8, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
plt.imshow(th3, cmap='gray')
<matplotlib.image.AxesImage at 0x208d0ee2070>
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
max_cntr_index = np.argmax(area_cntrs)
cnt = contours[max_cntr_index]
hull = cv2.convexHull(cnt)
if len(contours) < 2:
print("Errors")
cv2.drawContours(hull_img, [hull], 0, 125, 2)
plt.imshow(hull_img, cmap='gray')
<matplotlib.image.AxesImage at 0x208d0f92f10>
# hull_img = np.zeros_like(g_img_8, dtype=np.uint8)
# area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
# max_cntr_index = np.argmax(area_cntrs)
# cnt = contours[max_cntr_index]
# hull = cv2.convexHull(cnt)
# cv2.drawContours(hull_img, [hull], 0, 125, 2)
# plt.imshow(hull_img, cmap='gray')
ch_img_1 = np.zeros_like(img, dtype=np.uint8)
ch_img_1[:, :, 0] = hull_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_8
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
op[:10, :] = 0
op[-10:, :] = 0
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(op, cmap='gray')
plt.show()
# op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
# op[:10, :] = 0
# op[-10:, :] = 0
ch_img_1 = np.zeros_like(img, dtype=np.uint8)
ch_img_1[:, :, 0] = op
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_8+50
new_img = np.array((ch_img_1 + ch_img_2), dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
# ellipse fit below not terrific but not because of contour above is v. good
x, y = np.where(op > 1)
op_cnt = np.vstack([y, x]).T
ellipse = cv2.fitEllipse(op_cnt)
cv2.ellipse(img, ellipse, (2**15, 0, 0), 3)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(cv2.convertScaleAbs(img, alpha=(255/65535)) + 125)
plt.show()
ellipse
((2975.01220703125, 1009.5517578125), (5253.75927734375, 5340.55322265625), 97.18701934814453)
min_axis_len = 5253.75927734375
maj_axis_len = 5340.55322265625
ratio = min_axis_len/maj_axis_len
ratio
0.9837481358777039
def contour_points_v2(g_img):
kernel = np.ones((15, 15), np.uint8)
gradient = cv2.morphologyEx(g_img, cv2.MORPH_GRADIENT, kernel)
blur = cv2.GaussianBlur(gradient + 20, (5, 5), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
gradient_img = np.zeros_like(g_img, dtype=np.uint8)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 1:
return None, None
cv2.drawContours(gradient_img, contours, -1, 255, 2)
kernel = np.ones((25, 25),np.uint8)
gradient_img = cv2.dilate(gradient_img, kernel, iterations = 1)
hull_img = np.zeros_like(g_img, dtype=np.uint8)
blur = cv2.GaussianBlur(g_img, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 1:
return None, None
area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
max_cntr_index = np.argmax(area_cntrs)
cnt = contours[max_cntr_index]
hull = cv2.convexHull(cnt)
cv2.drawContours(hull_img, [hull], 0, 125, 2)
op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
# op[:10, :] = 0
# op[-10:, :] = 0
# kernel = np.ones((3, 3), np.uint8)
# opening = cv2.morphologyEx(op, cv2.MORPH_OPEN, kernel)
x, y = np.where(op > 1)
final_cnt = np.vstack([y, x]).T
# final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=int(final_cnt.shape[0]/10)))), axis=0)
# imgs = np.hstack([gradient_img, hull_img, opening])
if final_cnt.shape[0] > 4:
return op, final_cnt
else:
return None, None
# now preprocessing fixed for both types of images
# Ivo's implemetation of Ellipse class
class Ellipse_Fixed_Ratio:
C : np.array #homogeneous conic matrix x'*C*x=0
def __init__( self, C : np.array ):
self.C = C
@classmethod
def estimate_ellipsoid( cls, points : np.array ):
"""
C = estimate_ellipsoid_2D ( points )
----------------------------------
estimates a homogeneous ellipsoid C
x'*C*x = 0
from a number of point samples
(in a least squares sense)
format of points
[ x^1_1 . . . . . x^N_1 ]
[ . . ]
points = [ . . ]
[ . . ]
[ x^1_n . . . . . x^N_n ]
[ 1 . . . . . 1 ]
where each point is element R^n
and there are N points, the points
have to be in homogeneous coordinates.
"""
n = points.shape[0]
# n is 3
num_entries = int((n+1) * n/2) - 2 # fixing ellipse maj/min axis length ratio
if points.shape[1] < num_entries-1:
print('insufficient number of points, need at least {}'.format( num_entries-1 ) )
# -------------- estimate ellipse ---------------
# assemble linear system
A = np.zeros( [points.shape[ 1 ], num_entries ] )
#for all points
# print(points.shape[ 1 ])
for i in range( 0, points.shape[ 1 ] ):
#extract point
x = np.array([points[:, i]]).T
X = np.matmul( x, x.T )
X2 = X.T + X - np.diag( np.diag( X ) )
row = X2[ np.where( np.triu( np.ones( X2.shape ) ) ) ]
row[0] = row[0] + row[3]
row = np.delete(row, 3)
row = np.delete(row, 1)
A[ i, : ] = row
#solve for null space
U,D,V = np.linalg.svd( A )
# print(U, D, V)
ell = V[-1,:]
ell = np.insert(ell, 1, 0)
ell = np.insert(ell, 3, ell[0])
# ell[0] = 0.9837*ell[0]
print('Ell:', ell)
#create conic matrix
Cd = np.zeros( X2.shape )
Cd[ np.where( np.triu( np.ones( X2.shape ) ) ) ] = ell
C = Cd + Cd.T - np.diag( np.diag ( Cd ) )
print(f'C: {C}')
if Ellipse_Fixed_Ratio.is_hyperboloid( C ):
print('Warning: result is a hyperboloid!')
if Ellipse_Fixed_Ratio.is_paraboloid( C ):
print('Warning: result is a paraboloid!')
return C
@classmethod
def cov_mean_to_ellipsoid( cls, cov, mean, k ):
"""
C = cov_mean_to_ellipsoid_homog ( Cov, Mean, k )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from covariance matrix and mean value
of a gaussian distribution
k is a confidence factor determibed by the chi^2
distribution. The probability that x is inside
the ellipsoid is P_{chi^2}(k,2). Some values
for k are 2.41 (70%) and k=5.99 (95%).
cov : np.array - covariance matrix (non-homogeneous)
mean : np.array - mean value (vector)
k : float - confidence factor, see above
"""
#column vector, assume we get a 1D array
mean = np.array( [mean] ).T
SIGMAinv = np.linalg.inv( cov )
lastentry = np.matmul( mean.T, np.matmul( SIGMAinv, mean ) ) - k
C = np.vstack( [ np.hstack( [ SIGMAinv, -np.matmul(SIGMAinv, mean)] ),
np.hstack( [ -np.matmul(mean.T,SIGMAinv), lastentry] ) ] )
return C
@classmethod
def srt_to_ellipsoid( cls, S, R, T ):
"""
C = srt_to_ellipsoid_homog ( S, R, T )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from translation, rotation and scale parts as
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
see: Ivo Ihrke: "Some Notes on Ellipses", TechReport, 2004
"""
#function C = srt_to_ellipsoid_homog ( S, R, T )
invT = np.linalg.inv( T )
invR = R.T
invS = np.linalg.inv( S )
left = np.matmul( invT.T, np.matmul( invR.T, invS.T ) )
right = np.matmul( invS, np.matmul( invR, invT ) )
C = np.matmul( left, right )
return C
@classmethod
def decompose_ellipsoid( cls, C ):
"""
[S,R,T] = decompose_ellipsoid_homog ( C )
----------------------------------
decomposes a homogeneous ellipsoid C
x'*C*x = 0
into translation, rotation and scale parts
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
Note that C \neq (T^-1 * R^-1 * S^-1) ( S * R * T )
in some cases because the orientation of
the ellipsoid's major axis is not fixed.
The function generates a right handed system.
"""
#print(C)
C_ih = C[0:-1,0:-1]
#print(C_ih)
t = -np.matmul( np.linalg.pinv( C_ih ), C[0:-1,-1] )
T = np.eye( C.shape[ 0 ] )
T[0:-1, -1 ] = t
U,D,V = np.linalg.svd( C_ih )
R = np.eye( C.shape[ 0 ] )
R[0:-1,0:-1] = V
#S = R' * T' * C * T * R;
S = np.matmul( R.T, np.matmul( T.T, np.matmul( C, np.matmul( T, R ) ) ) )
S = S / S[-1,-1]
S = np.sqrt( np.diag(1/abs(np.diag(S))))
return S, R, T
@classmethod
def is_ellipsoid( cls, C ):
"""
bool = is_ellipsoid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is an ellipsoid
"""
#function val = conic_is_ellipsoid_homog ( C )
detC = np.linalg.det( C[0:-1,0:-1] )
return detC > 0
@classmethod
def is_hyperboloid( cls, C ):
"""
bool = is_hyperboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a hyperboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
return detC < 0
@classmethod
def is_paraboloid( cls, C ):
"""
bool = is_paraboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a paraboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
#maybe with EPS ?
return detC == 0
def get_center( self ):
S,R,T = self.decompose_ellipsoid( self.C )
return [ T[2,0], T[2,1] ]
def plot( self, ax, **kwargs ):
"""
plot_ellipsoid_2D ( C, style, N )
----------------------------------
plots a homogeneous ellipsoid C
x'*C*x = 0
style is a matlab plot style, e.g. 'b-' (default)
N is the number of polygon segments to
approximate the ellipse
"""
#function plot_ellipsoid_homog_2d ( C, style, N )
if not 'style' in locals():
style = 'b-'
if not 'N' in locals():
N = 100
angles = np.linspace( -np.pi, np.pi, N )
# generate a unit circle or a unit hyperbola
if ( Ellipse_Fixed_Ratio.is_ellipsoid( self.C ) ):
#uc = [ cos(x); sin(x); ones(size(x)) ];
x = np.cos( angles )
y = np.sin( angles )
w = np.ones( angles.shape )
uc = np.vstack([x,y,w])
else:
EPS = 0.1
x = np.linspace( -np.pi/2+EPS, np.pi/2-EPS, int( N/2 ) )#-(pi/2-EPS):((pi-2*EPS)/(N/2)):(pi/2-EPS);
x2 = np.linspace( np.pi/2+EPS, 3*np.pi/2-EPS, int( N/2 ) ) #x2 = (pi/2+EPS):((pi-2*EPS)/(N/2)):(3*pi/2-EPS);
uc = np.vstack( [ 1/np.cos(x), np.tan(x), np.ones( x.shape ) ] )
uc2 = np.vstack( [ 1/np.cos(x2), np.tan(x2), np.ones( x2.shape ) ] )
#decompose
S,R,T = Ellipse_Fixed_Ratio.decompose_ellipsoid( self.C )
#transform unit circle or parabola
transform = np.matmul( T, np.matmul( R, S ) )
ell = np.matmul( transform, uc )
if 'uc2' in locals():
ell2 = np.matmul( transform, uc2 )
ax.plot( ell2[0,:], ell2[1,:], style )
#plot conic
ax.plot( ell[0,:], ell[1,:], style )
#plot coordinate system
tx = T[0,2]
ty = T[1,2]
# print(tx, ty)
smin = S[0,0]
smaj = S[1,1]
ax.plot( np.array( [ tx, tx+smin*R[0,0] ] ), np.array( [ ty, ty+smin*R[1,0] ] ), 'r' )
ax.plot( np.array( [ tx, tx+smaj*R[0,1] ] ), np.array( [ ty, ty+smaj*R[1,1] ] ), 'b' )
ax.set_aspect('equal', adjustable='box')
return (tx, ty, smin, smaj)
def fit_ellipse_2(path_dir_w_time, ifHull=False, date_now=""):
""" Find and fit an ellipse in an image.
it should be run on a normalized image I (normalization should exclude dead pixels which are artificially maximum)
Parameters
----------
path_dir_w_time : str
_description_
date_now : str
_description_
Returns
-------
Ellipse
Ellipse in an image.
"""
I_temp = plt.imread(path_dir_w_time)[:,:,1].astype(float)/255
I = old_rplc_dead_px_med(I_temp)
ellipsoid_points=[]
if ifHull == False:
#every Nth scanline should be sampled
START_SCANLINE = 50 #there are some weird effects on the zero scanline that lead to detection of erroneous points
DELTA_SCANLINE = 100
#need to filter enough such that the big jump into and out of the disc has indeed the largest magnitude derivative
FILTERWIDTH=250
x=np.linspace(-4,4,FILTERWIDTH)
gauss=np.exp(-x**2)
plt.figure(figsize=[18,12])
for sl in range( START_SCANLINE, I.shape[0]-START_SCANLINE, DELTA_SCANLINE ):
#print(f'Checking Scanline {sl} ...' )
#line = self.I[sl]
line = I[sl, :]
filtered_line = sgn.convolve(line,gauss,mode='same')
#crude filter for line intersecting the ring structure
if np.amax(filtered_line) > 0.9:
#print("- scanline intersection ")
filtered_deriv = np.diff( filtered_line )
#crude extraction - this will probably need improvement
max_pos = np.argmax( filtered_deriv )
min_pos = np.argmin( filtered_deriv )
ellipsoid_points.append([max_pos,sl])
ellipsoid_points.append([min_pos,sl]) #add two points on perimeter
else:
#print("- no intersection ")
pass
print(ellipsoid_points)
else:
img = np.array(I*255, dtype=np.uint8)
cnt_img, final_cnt = contour_points_v2(img)
if final_cnt is None:
return None, None
elif final_cnt.shape[0] < 5:
return None, None
else:
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
# print(ellipsoid_points)
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
# print(ellipsoid_points)
#plt.clf()
#plt.imshow(self.I)
#plt.show()
for i in range( 0, ellipsoid_points.shape[1] ):
plt.plot(ellipsoid_points[0,i],ellipsoid_points[1,i],'rx')
E = Ellipse_Fixed_Ratio( Ellipse_Fixed_Ratio.estimate_ellipsoid( ellipsoid_points ) )
S,R,T=Ellipse_Fixed_Ratio.decompose_ellipsoid(E.C)
print(f'S:{S}, \nR: {R}, \nT: {T}')
plt.plot(T[0,2],T[1,2],'gx')
plt.imshow(I)
center = E.plot(plt.gca())
plt.title(f'ellipsoide fit {1}')
# plt.savefig(f"{path_dir_w_time}/{date_now}_ledpos_{self.x_pos}_{self.y_pos}_{self.channel}_ellipsoide_fit.png")
return E, center
def old_rplc_dead_px_med(I_temp):
"""Replace dead pixels by median
Parameters
----------
I_temp : ndarray
image to be processed
Returns
-------
ndarray
clean image where dead pixels are replaced by median
"""
Imed = sgn.medfilt2d( I_temp, kernel_size=5 )
maxI = np.amax( Imed ) #normalize robustly
I_temp = I_temp / maxI
#replace dead pixels by median filtered values
I_temp[ I_temp > 1 ] = Imed[ I_temp > 1 ]
return I_temp
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"
# dir_path = "Z:/CSE/CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_08_25/2023_08_25_09_31_25/"
# file_path = "2023_08_25_09_31_25_img_iso_100_shutter_2_x_15_y_15_r_0_g_1_b_0.tiff"
E, center = fit_ellipse(dir_path + file_path, ifHull=True)
plt.show()
print(center)
[[1.440e+02 1.450e+02 1.460e+02 1.440e+02 1.450e+02 1.460e+02 1.440e+02 1.450e+02 1.460e+02 4.557e+03 4.558e+03 4.559e+03 4.556e+03 4.557e+03 4.558e+03 4.559e+03 4.556e+03 4.557e+03 4.558e+03 4.559e+03 4.556e+03 4.557e+03 4.558e+03 4.480e+03 4.481e+03 4.482e+03 4.480e+03 4.481e+03 4.482e+03 4.480e+03 4.481e+03 4.482e+03 4.325e+03 4.326e+03 4.327e+03 4.325e+03 4.326e+03 4.327e+03 4.325e+03 4.326e+03 4.327e+03 3.509e+03 3.510e+03 3.511e+03 3.501e+03 3.502e+03 3.503e+03 3.504e+03 3.505e+03 3.509e+03 3.510e+03 3.511e+03 3.496e+03 3.497e+03 3.498e+03 3.499e+03 3.500e+03 3.501e+03 3.502e+03 3.503e+03 3.504e+03 3.505e+03 3.509e+03 3.510e+03 3.511e+03 3.492e+03 3.493e+03 3.494e+03 3.495e+03 3.496e+03 3.497e+03 3.498e+03 3.499e+03 3.500e+03 3.501e+03 3.502e+03 3.503e+03 3.504e+03 3.505e+03 3.486e+03 3.487e+03 3.488e+03 3.492e+03 3.493e+03 3.494e+03 3.495e+03 3.496e+03 3.497e+03 3.498e+03 3.499e+03 3.500e+03 3.486e+03 3.487e+03 3.488e+03 3.492e+03 3.493e+03 3.494e+03 3.495e+03 3.486e+03 3.487e+03 3.488e+03] [1.791e+03 1.791e+03 1.791e+03 1.792e+03 1.792e+03 1.792e+03 1.793e+03 1.793e+03 1.793e+03 3.192e+03 3.192e+03 3.192e+03 3.193e+03 3.193e+03 3.193e+03 3.193e+03 3.194e+03 3.194e+03 3.194e+03 3.194e+03 3.195e+03 3.195e+03 3.195e+03 3.265e+03 3.265e+03 3.265e+03 3.266e+03 3.266e+03 3.266e+03 3.267e+03 3.267e+03 3.267e+03 3.367e+03 3.367e+03 3.367e+03 3.368e+03 3.368e+03 3.368e+03 3.369e+03 3.369e+03 3.369e+03 3.730e+03 3.730e+03 3.730e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.736e+03 3.736e+03 3.736e+03] [1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00]] 3 Warning: result is a hyperboloid!
(5272.143308579376, 3112.6301639181233)
E, center = fit_ellipse_2(dir_path + file_path, ifHull=True)
plt.show()
print(center)
Ell: [ 5.23202987e-07 0.00000000e+00 -1.44310664e-03 5.14674778e-07 -6.28543710e-04 9.99998761e-01] C: [[ 5.23202987e-07 0.00000000e+00 -1.44310664e-03] [ 0.00000000e+00 5.14674778e-07 -6.28543710e-04] [-1.44310664e-03 -6.28543710e-04 9.99998761e-01]] S:[[2.67648633e+03 0.00000000e+00 0.00000000e+00] [0.00000000e+00 2.69857004e+03 0.00000000e+00] [0.00000000e+00 0.00000000e+00 1.00000000e+00]], R: [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]], T: [[1.00000000e+00 0.00000000e+00 2.75821559e+03] [0.00000000e+00 1.00000000e+00 1.22124444e+03] [0.00000000e+00 0.00000000e+00 1.00000000e+00]]
(2758.2155895125884, 1221.2444375066455, 2676.486329255531, 2698.5700358857857)
center[2]/center[3]
0.9918165152889924
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, 195, axis=1))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(20, 20))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, 0), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, 0), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_19_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, 0), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -195), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -195), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -210), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_19_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -200), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (177, -180), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (187, -205), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, -200), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
g_img_2 = g_img.copy()
ellipse = ((2755+390+15, 1210-390),
(5322, 5200),
0)
# 17,17 -> not a good fit: more elliptical
# ((2755+390+15, 1210-390),
# (5322, 5200),
# 0)
# 13,13
# ((2755-390, 1210+390),
# (5322, 5322),
# 0)
# 15, 17
# ((2755, 1210-380),
# (5322, 5200),
# 0)
# 15, 13
# ((2755, 1210+400/390),
# (5320, 5320),
# 0)
# 13,15:
# ((2765-410, 1210-15),
# (5320, 5320),
# 0)
# 15,15:
# ((2765, 1210),
# (5320, 5320),
# 0)
# 17,15:
# ((2765+390, 1210-20),
# (5320, 5320),
# 0)
cv2.ellipse(g_img_2, ellipse, 255, 4)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(g_img_2, cmap='gray')
plt.show()
def contour_points_v2(g_img):
kernel = np.ones((15, 15), np.uint8)
gradient = cv2.morphologyEx(g_img, cv2.MORPH_GRADIENT, kernel)
blur = cv2.GaussianBlur(gradient + 20, (5, 5), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
gradient_img = np.zeros_like(g_img, dtype=np.uint8)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 1:
return None
cv2.drawContours(gradient_img, contours, -1, 255, 2)
kernel = np.ones((25, 25),np.uint8)
gradient_img = cv2.dilate(gradient_img, kernel, iterations = 1)
hull_img = np.zeros_like(g_img, dtype=np.uint8)
blur = cv2.GaussianBlur(g_img, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 1:
return None
area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
max_cntr_index = np.argmax(area_cntrs)
cnt = contours[max_cntr_index]
hull = cv2.convexHull(cnt)
cv2.drawContours(hull_img, [hull], 0, 125, 2)
op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
# op[:10, :] = 0
# op[-10:, :] = 0
# kernel = np.ones((3, 3), np.uint8)
# opening = cv2.morphologyEx(op, cv2.MORPH_OPEN, kernel)
# final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=int(final_cnt.shape[0]/10)))), axis=0)
# imgs = np.hstack([gradient_img, hull_img, opening])
if final_cnt.shape[0] > 4:
return op
else:
return None
# Ivo's implemetation of Ellipse class
class Ellipse:
C : np.array #homogeneous conic matrix x'*C*x=0
def __init__( self, C : np.array ):
self.C = C
@classmethod
def estimate_ellipsoid( cls, points : np.array ):
"""
C = estimate_ellipsoid_2D ( points )
----------------------------------
estimates a homogeneous ellipsoid C
x'*C*x = 0
from a number of point samples
(in a least squares sense)
format of points
[ x^1_1 . . . . . x^N_1 ]
[ . . ]
points = [ . . ]
[ . . ]
[ x^1_n . . . . . x^N_n ]
[ 1 . . . . . 1 ]
where each point is element R^n
and there are N points, the points
have to be in homogeneous coordinates.
"""
n = points.shape[0]
num_entries = int((n+1) * n/2)
if points.shape[1] < num_entries-1:
print('insufficient number of points, need at least {}'.format( num_entries-1 ) )
# -------------- estimate ellipse ---------------
# assemble linear system
A = np.zeros( [points.shape[ 1 ], num_entries ] )
#for all points
for i in range( 0, points.shape[ 1 ] ):
#extract point
x = np.array([points[:, i]]).T
#print(x)
X = np.matmul( x, x.T )
#print(X)
X2 = X.T + X - np.diag( np.diag( X ) )
#print(X2)
row = X2[ np.where( np.triu( np.ones( X2.shape ) ) ) ]
#print(row)
A[ i, : ] = row
#print(A)
#solve for null space
U,D,V = np.linalg.svd( A )
ell = V[-1,:]
#create conic matrix
Cd = np.zeros( X2.shape )
Cd[ np.where( np.triu( np.ones( X2.shape ) ) ) ] = ell
C = Cd + Cd.T - np.diag( np.diag ( Cd ) )
#detC = np.linalg.det( C[0:-1,0:-1] )
#print(detC)
if Ellipse.is_hyperboloid( C ):
print('Warning: result is a hyperboloid!')
if Ellipse.is_paraboloid( C ):
print('Warning: result is a paraboloid!')
return C
@classmethod
def cov_mean_to_ellipsoid( cls, cov, mean, k ):
"""
C = cov_mean_to_ellipsoid_homog ( Cov, Mean, k )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from covariance matrix and mean value
of a gaussian distribution
k is a confidence factor determibed by the chi^2
distribution. The probability that x is inside
the ellipsoid is P_{chi^2}(k,2). Some values
for k are 2.41 (70%) and k=5.99 (95%).
cov : np.array - covariance matrix (non-homogeneous)
mean : np.array - mean value (vector)
k : float - confidence factor, see above
"""
#column vector, assume we get a 1D array
mean = np.array( [mean] ).T
SIGMAinv = np.linalg.inv( cov )
lastentry = np.matmul( mean.T, np.matmul( SIGMAinv, mean ) ) - k
C = np.vstack( [ np.hstack( [ SIGMAinv, -np.matmul(SIGMAinv, mean)] ),
np.hstack( [ -np.matmul(mean.T,SIGMAinv), lastentry] ) ] )
return C
@classmethod
def srt_to_ellipsoid( cls, S, R, T ):
"""
C = srt_to_ellipsoid_homog ( S, R, T )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from translation, rotation and scale parts as
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
see: Ivo Ihrke: "Some Notes on Ellipses", TechReport, 2004
"""
#function C = srt_to_ellipsoid_homog ( S, R, T )
invT = np.linalg.inv( T )
invR = R.T
invS = np.linalg.inv( S )
left = np.matmul( invT.T, np.matmul( invR.T, invS.T ) )
right = np.matmul( invS, np.matmul( invR, invT ) )
C = np.matmul( left, right )
return C
@classmethod
def decompose_ellipsoid( cls, C ):
"""
[S,R,T] = decompose_ellipsoid_homog ( C )
----------------------------------
decomposes a homogeneous ellipsoid C
x'*C*x = 0
into translation, rotation and scale parts
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
Note that C \neq (T^-1 * R^-1 * S^-1) ( S * R * T )
in some cases because the orientation of
the ellipsoid's major axis is not fixed.
The function generates a right handed system.
"""
#print(C)
C_ih = C[0:-1,0:-1]
#print(C_ih)
t = -np.matmul( np.linalg.pinv( C_ih ), C[0:-1,-1] )
T = np.eye( C.shape[ 0 ] )
T[0:-1, -1 ] = t
U,D,V = np.linalg.svd( C_ih )
R = np.eye( C.shape[ 0 ] )
R[0:-1,0:-1] = V
#S = R' * T' * C * T * R;
S = np.matmul( R.T, np.matmul( T.T, np.matmul( C, np.matmul( T, R ) ) ) )
S = S / S[-1,-1]
S = np.sqrt( np.diag(1/abs(np.diag(S))))
return S, R, T
@classmethod
def is_ellipsoid( cls, C ):
"""
bool = is_ellipsoid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is an ellipsoid
"""
#function val = conic_is_ellipsoid_homog ( C )
detC = np.linalg.det( C[0:-1,0:-1] )
return detC > 0
@classmethod
def is_hyperboloid( cls, C ):
"""
bool = is_hyperboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a hyperboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
return detC < 0
@classmethod
def is_paraboloid( cls, C ):
"""
bool = is_paraboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a paraboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
#maybe with EPS ?
return detC == 0
def get_center( self ):
S,R,T = self.decompose_ellipsoid( self.C )
return [ T[2,0], T[2,1] ]
def plot( self, ax, **kwargs ):
"""
plot_ellipsoid_2D ( C, style, N )
----------------------------------
plots a homogeneous ellipsoid C
x'*C*x = 0
style is a matlab plot style, e.g. 'b-' (default)
N is the number of polygon segments to
approximate the ellipse
"""
#function plot_ellipsoid_homog_2d ( C, style, N )
if not 'style' in locals():
style = 'b-'
if not 'N' in locals():
N = 100
angles = np.linspace( -np.pi, np.pi, N )
# generate a unit circle or a unit hyperbola
if ( Ellipse.is_ellipsoid( self.C ) ):
#uc = [ cos(x); sin(x); ones(size(x)) ];
x = np.cos( angles )
y = np.sin( angles )
w = np.ones( angles.shape )
uc = np.vstack([x,y,w])
else:
EPS = 0.1
x = np.linspace( -np.pi/2+EPS, np.pi/2-EPS, int( N/2 ) )#-(pi/2-EPS):((pi-2*EPS)/(N/2)):(pi/2-EPS);
x2 = np.linspace( np.pi/2+EPS, 3*np.pi/2-EPS, int( N/2 ) ) #x2 = (pi/2+EPS):((pi-2*EPS)/(N/2)):(3*pi/2-EPS);
uc = np.vstack( [ 1/np.cos(x), np.tan(x), np.ones( x.shape ) ] )
uc2 = np.vstack( [ 1/np.cos(x2), np.tan(x2), np.ones( x2.shape ) ] )
#decompose
S,R,T = Ellipse.decompose_ellipsoid( self.C )
#transform unit circle or parabola
transform = np.matmul( T, np.matmul( R, S ) )
ell = np.matmul( transform, uc )
if 'uc2' in locals():
ell2 = np.matmul( transform, uc2 )
ax.plot( ell2[0,:], ell2[1,:], style )
#plot conic
ax.plot( ell[0,:], ell[1,:], style )
#plot coordinate system
tx = T[0,2]
ty = T[1,2]
# print(tx, ty)
smin = S[0,0]
smaj = S[1,1]
ax.plot( np.array( [ tx, tx+smin*R[0,0] ] ), np.array( [ ty, ty+smin*R[1,0] ] ), 'r' )
ax.plot( np.array( [ tx, tx+smaj*R[0,1] ] ), np.array( [ ty, ty+smaj*R[1,1] ] ), 'b' )
ax.set_aspect('equal', adjustable='box')
return (tx, ty)
def fit_ellipse_3(img_name, img_8, g_img_8, ifHull=False, date_now="", method="ivo", save_path=None):
""" Find and fit an ellipse in an image.
it should be run on a normalized image I (normalization should exclude dead pixels which are artificially maximum)
Parameters
----------
path_dir_w_time : str
_description_
date_now : str
_description_
Returns
-------
Ellipse
Ellipse in an image.
"""
ellipsoid_points=[]
img = np.array(g_img_8, dtype=np.uint8)
op = contour_points_v2(img)
op[:20, :] = 0
op[-20:, :] = 0
op[:, :20] = 0
op[:, -20:] = 0
# plt.imshow(op)
# plt.show()
x, y = np.where(op > 1)
final_cnt = np.vstack([y, x]).T
if final_cnt.shape[0] > 10000:
final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=10000))), axis=0)
if final_cnt is None:
return None, None, None
else:
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
# print(ellipsoid_points)
is_det_ellipse = True
if method == "ivo":
E = Ellipse( Ellipse.estimate_ellipsoid( ellipsoid_points ) )
if Ellipse.is_ellipsoid( E.C ):
S,R,T = Ellipse.decompose_ellipsoid(E.C)
tx = T[0,2]
ty = T[1,2]
smin = S[0,0]
smaj = S[1,1]
# f = E.C[-1, -1]
f = np.arctan2(R[1,0], R[0, 0])
# print(theta)
# print(E.C)
else:
is_det_ellipse = False
else:
(tx, ty), (smin, smaj), f = cv2.fitEllipse(ellipsoid_points[:2, :].T.astype(np.int32))
smin, smaj = smin/2, smaj/2
# print((tx, ty), (smin, smaj), f)
#
if is_det_ellipse:
for i in range( 0, ellipsoid_points.shape[1] ):
draw_cross(img_8, (ellipsoid_points[0,i],ellipsoid_points[1,i]), (0, 0, 255), length=20, thickness=5)
cv2.ellipse(img_8, (int(tx), int(ty)), (int(smin), int(smaj)), f, 0, 360, (255, 0, 0), 5)
draw_cross(img_8, (tx, ty), (0, 255, 0))
# plt.imshow(img_8)
# plt.show()
cv2.imshow('Resized Image', cv2.resize(img_8, (1080, 720), interpolation = cv2.INTER_LINEAR))
cv2.waitKey(0)
cv2.destroyAllWindows()
if save_path is not None:
cv2.imwrite(f'{save_path}/{img_name[:-5]}.png', img_8)
return (int(tx), int(ty)), (int(smin), int(smaj)), f
else:
return None, None, None
cv2.fitEllipse(ellipsoid_points[:2, :].T.astype(np.int32))
((219.3291473388672, 104.70303344726562), (513.478759765625, 532.2569580078125), 53.40561294555664)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
center, axis, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True, method='ivo')
[[ 7.23449895e-07 7.32307094e-09 -1.69863270e-03] [ 7.32307094e-09 7.46098008e-07 -4.49513750e-04] [-1.69863270e-03 -4.49513750e-04 -9.99998456e-01]]
def draw_cross(img, cord, color, length=50, thickness=10):
x, y = int(cord[0]), int(cord[1])
cv2.line(img, (x-length, y-length), (x+length, y+length), color, thickness)
cv2.line(img, (x+length, y-length), (x-length, y+length), color, thickness)
return img
# # grid based on thumbnail images
# fit_range = np.zeros([7, 7], dtype=np.uint8)
# center_cords_ivo = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_ivo = np.zeros([7, 7, 2], dtype=np.float64)
# for x in range(12, 19):
# for y in range(12, 19):
# dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
# file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
# print('x:', x, 'y:', y, file_path)
# img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
# _, g_img, _ = cv2.split(img)
# img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
# g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
# center, axis_size, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True)
# if center is not None:
# good_fit = input("Is it a good fit?")
# fit_range[x-12, y-12] = int(good_fit)
# center_cords_ivo[x-12, y-12] = center
# axis_lens_ivo[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff Is it a good fit?5 x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff Is it a good fit?3 x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff Is it a good fit?4 x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff Is it a good fit?2 x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff Is it a good fit?2 x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff Warning: result is a hyperboloid! x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff Is it a good fit?3 x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff Is it a good fit?3 x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff Warning: result is a hyperboloid! x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff Is it a good fit?2 x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff Is it a good fit?3 x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff Is it a good fit?3 x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff Warning: result is a paraboloid! x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff Is it a good fit?2 x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff Warning: result is a paraboloid! x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff Is it a good fit?3 x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff Is it a good fit?2 x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff Is it a good fit?3 x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff Is it a good fit?3 x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff Warning: result is a paraboloid! x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff Is it a good fit?3 x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff Is it a good fit?3 x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff Is it a good fit?5 x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff Warning: result is a paraboloid! x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff Is it a good fit?3 x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff Is it a good fit?2 x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff Is it a good fit?5
plt.figure(figsize=(15, 10))
center_leds = center_cords_ivo[3, 3]
y_unit_len = center_cords_ivo[3, 2][1] - center_cords_ivo[3, 3][1]
x_unit_len = center_cords_ivo[4, 3][0] - center_cords_ivo[3, 3][0]
center_leds, y_unit_len, x_unit_len
x_axis = [center_leds[0] - i*x_unit_len for i in range(-6, 6)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-5, 5)]
xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')
x_cords = [x if x != 0.0 else np.nan for x in center_cords_ivo.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_ivo.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range.reshape(-1)]
plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()
count = 0
for x in range(12, 19):
for y in range(12, 19):
plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
count += 1
plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
axis_lens_pro = axis_lens_ivo.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')
xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])
cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')
cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')
plt.show()
axis_lens_pro = axis_lens_ivo.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
# # grid based on thumbnail images
# fit_range_cv2 = np.zeros([7, 7], dtype=np.uint8)
# center_cords_cv2 = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_cv2 = np.zeros([7, 7, 2], dtype=np.float64)
# save_dir = "C:/Users/Kazim/Desktop/Overall/VignettingCircleDetection/310823_ellipse_cv2_fit_tiff_size"
# for x in range(12, 19):
# for y in range(12, 19):
# dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
# file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
# print('x:', x, 'y:', y, file_path)
# img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
# _, g_img, _ = cv2.split(img)
# img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
# g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
# center, axis_size, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True, method='opencv', save_path=save_dir)
# if center is not None:
# good_fit = input("Is it a good fit?")
# fit_range_cv2[x-12, y-12] = int(good_fit)
# center_cords_cv2[x-12, y-12] = center
# axis_lens_cv2[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff Is it a good fit?5 x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff Is it a good fit?2 x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff Is it a good fit?5 x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff Is it a good fit?2 x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff Is it a good fit?5
plt.figure(figsize=(15, 10))
center_leds = center_cords_cv2[3, 3]
y_unit_len = center_cords_cv2[3, 2][1] - center_cords_cv2[3, 3][1]
x_unit_len = center_cords_cv2[4, 3][0] - center_cords_cv2[3, 3][0]
center_leds, y_unit_len, x_unit_len
x_axis = [center_leds[0] - i*x_unit_len for i in range(-4, 4)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-4, 4)]
xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')
x_cords = [x if x != 0.0 else np.nan for x in center_cords_cv2.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_cv2.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range_cv2.reshape(-1)]
plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()
count = 0
for x in range(12, 19):
for y in range(12, 19):
plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
count += 1
plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
axis_lens_pro = axis_lens_cv2.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')
xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])
cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')
cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')
plt.show()
axis_lens_pro = axis_lens_cv2.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
# fixing A, B, C, F
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
center, axis, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True, method='ivo') # good fit
[[-3.24326055e-07 4.75306325e-10 9.59646019e-04] [ 4.75306325e-10 -3.27813554e-07 3.93645616e-04] [ 9.59646019e-04 3.93645616e-04 -9.99999462e-01]]
np.cos()
img = np.array(g_img_8, dtype=np.uint8)
op = contour_points_v2(img)
op[:20, :] = 0
op[-20:, :] = 0
op[:, :20] = 0
op[:, -20:] = 0
x, y = np.where(op > 1)
final_cnt = np.vstack([y, x]).T
if final_cnt.shape[0] > 10000:
final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=10000))), axis=0)
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points
array([[ 529., 530., 531., ..., 2178., 2179., 2180.],
[ 93., 93., 93., ..., 3732., 3732., 3732.]])
b = -3.24326055e-07*x**2 + -3.27813554e-07*y**2 + f
A = np.array([x, y]).T
A, -b
(array([[ 93, 529],
[ 93, 530],
[ 93, 531],
...,
[3732, 2178],
[3732, 2179],
[3732, 2180]], dtype=int64),
array([1.09454023, 1.09488739, 1.0952352 , ..., 7.07219959, 7.07362787,
7.07505681]))
np.linalg.lstsq(A, -b, rcond=None)
(array([0.00077151, 0.00194168]), array([0.65526845]), 2, array([338927.99799414, 121566.05420793]))
9.59646019e-04*2, 3.93645616e-04*2
(0.001919292038, 0.000787291232)
9.26610761e-04*2, -9.39963645e-04*2
(0.001853221522, -0.00187992729)
b = 4.75306325e-10*x*y + -9.99999462e-01
A = np.array([x**2, y**2, x, y]).T
A, -b
(array([[ 8649, 279841, 93, 529],
[ 8649, 280900, 93, 530],
[ 8649, 281961, 93, 531],
...,
[13927824, 4743684, 3732, 2178],
[13927824, 4748041, 3732, 2179],
[13927824, 4752400, 3732, 2180]], dtype=int64),
array([0.99997608, 0.99997603, 0.99997599, ..., 0.99613603, 0.99613426,
0.99613248]))
A, B, D, E = np.linalg.lstsq(A, -b, rcond=None)[0]
A, B, D/2, E/2
(-3.268969844369404e-07, -3.236289044146954e-07, 0.0003930310412070922, 0.0009580904040238742)
def fit_ellipse_4(img_name, img_8, g_img_8, save_path=None):
""" Find ellipse center position vector where A=-3.24326055e-07, B=0, C=-3.27813554e-07, F=-9.99999462e-01
Parameters
----------
Returns
-------
Ellipse
Ellipse in an image.
"""
ellipsoid_points=[]
img = np.array(g_img_8, dtype=np.uint8)
op = contour_points_v2(img)
op[:20, :] = 0
op[-20:, :] = 0
op[:, :20] = 0
op[:, -20:] = 0
x, y = np.where(op > 0)
final_cnt = np.vstack([y, x]).T
if final_cnt is None:
return None, None, None
else:
if final_cnt.shape[0] > 10000:
final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=10000))), axis=0)
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
is_det_ellipse = True
# C = -3.24326055e-07
B = 4.75306325e-10
# A = -3.27813554e-07
F = -9.99999462e-01*2
# b = A*x**2 + B*x*y + C*y**2 + F
b = B*x*y + F
# eqs = np.array([x**2, y**2, x, y]).T
# C, A, E, D = np.linalg.lstsq(eqs, -b, rcond=None)[0]/2
# matC = np.array([[A, B, D/2], [B, C, E/2], [D/2, E/2, F]])
eqs = np.array([x**2 + y**2, x, y]).T
A, E, D = np.linalg.lstsq(eqs, -b, rcond=None)[0]/2
matC = np.array([[A, B, D/2], [B, A, E/2], [D/2, E/2, F/2]])
E = Ellipse( matC )
# print(Ellipse.is_ellipsoid( E.C ))
print(matC)
# print(np.linalg.lstsq(eqs, -b, rcond=None))
if Ellipse.is_ellipsoid( E.C ):
S,R,T = Ellipse.decompose_ellipsoid(E.C)
tx = T[0,2]
ty = T[1,2]
smin = S[0,0]
smaj = S[1,1]
theta = np.arctan2(R[1,0], R[0, 0])
# print(theta)
# f = E.C[-1, -1]
# print(E.C)
else:
is_det_ellipse = False
if is_det_ellipse:
for i in range( 0, ellipsoid_points.shape[1] ):
draw_cross(img_8, (ellipsoid_points[0,i],ellipsoid_points[1,i]), (0, 0, 255), length=20, thickness=5)
cv2.ellipse(img_8, (int(tx), int(ty)), (int(smin), int(smaj)), theta, 0, 360, (255, 0, 0), 5)
draw_cross(img_8, (tx, ty), (0, 255, 0))
# plt.imshow(img_8)
# plt.show()
cv2.imshow('Resized Image', cv2.resize(img_8, (1080, 720), interpolation = cv2.INTER_LINEAR))
cv2.waitKey(0)
cv2.destroyAllWindows()
if save_path is not None:
cv2.imwrite(f'{save_path}/{img_name[:-5]}.png', img_8)
return (int(tx), int(ty)), (int(smin), int(smaj)), theta
else:
return None, None, None
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
center, axis_size, f = fit_ellipse_4(file_path, img_8, g_img_8, save_path=None)
[[-9.96262090e-07 4.75306325e-10 2.74601838e-03] [ 4.75306325e-10 -9.96262090e-07 7.62551767e-04] [ 2.74601838e-03 7.62551767e-04 -9.99999462e-01]]
# # grid based on thumbnail images
# fit_range_circle = np.zeros([7, 7], dtype=np.uint8)
# center_cords_circle = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_circle = np.zeros([7, 7, 2], dtype=np.float64)
# save_dir = "C:/Users/Kazim/Desktop/Overall/VignettingCircleDetection/310823_circle_fit_tiff_size"
# for x in range(12, 19):
# for y in range(12, 19):
# dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
# file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
# print('x:', x, 'y:', y, file_path)
# img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
# _, g_img, _ = cv2.split(img)
# img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
# g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
# center, axis_size, f = fit_ellipse_4(file_path, img_8, g_img_8, save_path=save_dir)
# if center is not None:
# good_fit = input("Is it a good fit?")
# fit_range_circle[x-12, y-12] = int(good_fit)
# center_cords_circle[x-12, y-12] = center
# axis_lens_circle[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff [[-1.51944827e-07 4.75306325e-10 4.28640374e-04] [ 4.75306325e-10 -1.51944827e-07 2.64667160e-04] [ 4.28640374e-04 2.64667160e-04 -9.99999462e-01]] Is it a good fit?5 x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff [[-8.05421194e-07 4.75306325e-10 1.84482464e-03] [ 4.75306325e-10 -8.05421194e-07 1.29403723e-03] [ 1.84482464e-03 1.29403723e-03 -9.99999462e-01]] Is it a good fit?3 x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff [[ 1.82876905e-06 4.75306325e-10 -3.95544562e-03] [ 4.75306325e-10 1.82876905e-06 -2.53086219e-03] [-3.95544562e-03 -2.53086219e-03 -9.99999462e-01]] Is it a good fit?1 x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff [[ 8.30978298e-07 4.75306325e-10 -1.78516376e-03] [ 4.75306325e-10 8.30978298e-07 -9.72636790e-04] [-1.78516376e-03 -9.72636790e-04 -9.99999462e-01]] Is it a good fit?2 x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff [[ 5.53879588e-07 4.75306325e-10 -1.18126017e-03] [ 4.75306325e-10 5.53879588e-07 -5.30531858e-04] [-1.18126017e-03 -5.30531858e-04 -9.99999462e-01]] Is it a good fit?2 x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff [[ 4.35008066e-07 4.75306325e-10 -9.22324007e-04] [ 4.75306325e-10 4.35008066e-07 -3.19552659e-04] [-9.22324007e-04 -3.19552659e-04 -9.99999462e-01]] Is it a good fit?2 x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff [[ 3.82984548e-07 4.75306325e-10 -8.12447687e-04] [ 4.75306325e-10 3.82984548e-07 -1.98597659e-04] [-8.12447687e-04 -1.98597659e-04 -9.99999462e-01]] Is it a good fit?2 x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff [[-5.58856538e-07 4.75306325e-10 1.33007273e-03] [ 4.75306325e-10 -5.58856538e-07 1.00255279e-03] [ 1.33007273e-03 1.00255279e-03 -9.99999462e-01]] Is it a good fit?1 x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff [[-9.23869595e-07 4.75306325e-10 2.18874479e-03] [ 4.75306325e-10 -9.23869595e-07 1.48275359e-03] [ 2.18874479e-03 1.48275359e-03 -9.99999462e-01]] Is it a good fit?1 x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff [[-2.27036215e-06 4.75306325e-10 5.37034775e-03] [ 4.75306325e-10 -2.27036215e-06 3.18468554e-03] [ 5.37034775e-03 3.18468554e-03 -9.99999462e-01]] Is it a good fit?1 x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff [[ 3.06726902e-06 4.75306325e-10 -7.19762759e-03] [ 4.75306325e-10 3.06726902e-06 -3.60634197e-03] [-7.19762759e-03 -3.60634197e-03 -9.99999462e-01]] Is it a good fit?2 x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff [[ 1.10581527e-06 4.75306325e-10 -2.58177670e-03] [ 4.75306325e-10 1.10581527e-06 -1.05802929e-03] [-2.58177670e-03 -1.05802929e-03 -9.99999462e-01]] Is it a good fit?2 x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff [[ 7.32686924e-07 4.75306325e-10 -1.70581000e-03] [ 4.75306325e-10 7.32686924e-07 -5.39313733e-04] [-1.70581000e-03 -5.39313733e-04 -9.99999462e-01]] Is it a good fit?2 x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff [[ 5.85509650e-07 4.75306325e-10 -1.36058178e-03] [ 4.75306325e-10 5.85509650e-07 -3.02023162e-04] [-1.36058178e-03 -3.02023162e-04 -9.99999462e-01]] Is it a good fit?2 x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff [[-3.66265272e-07 4.75306325e-10 9.41946218e-04] [ 4.75306325e-10 -3.66265272e-07 6.56238900e-04] [ 9.41946218e-04 6.56238900e-04 -9.99999462e-01]] Is it a good fit?1 x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff [[-4.88584245e-07 4.75306325e-10 1.25433906e-03] [ 4.75306325e-10 -4.88584245e-07 7.83008428e-04] [ 1.25433906e-03 7.83008428e-04 -9.99999462e-01]] Is it a good fit?1 x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff [[-7.02470516e-07 4.75306325e-10 1.80119649e-03] [ 4.75306325e-10 -7.02470516e-07 9.88621422e-04] [ 1.80119649e-03 9.88621422e-04 -9.99999462e-01]] Is it a good fit?1 x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff [[-1.26810745e-06 4.75306325e-10 3.23967294e-03] [ 4.75306325e-10 -1.26810745e-06 1.51545546e-03] [ 3.23967294e-03 1.51545546e-03 -9.99999462e-01]] Is it a good fit?2 x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff [[-3.00324635e-06 4.75306325e-10 7.67149042e-03] [ 4.75306325e-10 -3.00324635e-06 2.95940159e-03] [ 7.67149042e-03 2.95940159e-03 -9.99999462e-01]] Is it a good fit?2 x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff [[ 2.78928085e-06 4.75306325e-10 -7.07337439e-03] [ 4.75306325e-10 2.78928085e-06 -2.04190106e-03] [-7.07337439e-03 -2.04190106e-03 -9.99999462e-01]] Is it a good fit?2 x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff [[ 1.60175355e-06 4.75306325e-10 -4.06730316e-03] [ 4.75306325e-10 1.60175355e-06 -8.32599901e-04] [-4.06730316e-03 -8.32599901e-04 -9.99999462e-01]] Is it a good fit?2 x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff [[-2.66292824e-07 4.75306325e-10 7.36336097e-04] [ 4.75306325e-10 -2.66292824e-07 4.76573016e-04] [ 7.36336097e-04 4.76573016e-04 -9.99999462e-01]] Is it a good fit?1 x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff [[-3.25284852e-07 4.75306325e-10 8.98235317e-04] [ 4.75306325e-10 -3.25284852e-07 5.20749145e-04] [ 8.98235317e-04 5.20749145e-04 -9.99999462e-01]] Is it a good fit?1 x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff [[-4.03746681e-07 4.75306325e-10 1.11518083e-03] [ 4.75306325e-10 -4.03746681e-07 5.67327951e-04] [ 1.11518083e-03 5.67327951e-04 -9.99999462e-01]] Is it a good fit?1 x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff [[-5.21047365e-07 4.75306325e-10 1.43702280e-03] [ 4.75306325e-10 -5.21047365e-07 6.26811405e-04] [ 1.43702280e-03 6.26811405e-04 -9.99999462e-01]] Is it a good fit?1 x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff [[-7.03250637e-07 4.75306325e-10 1.93891868e-03] [ 4.75306325e-10 -7.03250637e-07 6.94528824e-04] [ 1.93891868e-03 6.94528824e-04 -9.99999462e-01]] Is it a good fit?2 x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff [[-9.96262090e-07 4.75306325e-10 2.74601838e-03] [ 4.75306325e-10 -9.96262090e-07 7.62551767e-04] [ 2.74601838e-03 7.62551767e-04 -9.99999462e-01]] Is it a good fit?2 x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff [[-1.36544789e-06 4.75306325e-10 3.75679371e-03] [ 4.75306325e-10 -1.36544789e-06 7.67007660e-04] [ 3.75679371e-03 7.67007660e-04 -9.99999462e-01]] Is it a good fit?2 x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff [[-2.06004470e-07 4.75306325e-10 6.09752267e-04] [ 4.75306325e-10 -2.06004470e-07 3.67305162e-04] [ 6.09752267e-04 3.67305162e-04 -9.99999462e-01]] Is it a good fit?1 x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff [[-2.38580862e-07 4.75306325e-10 7.06047598e-04] [ 4.75306325e-10 -2.38580862e-07 3.80434787e-04] [ 7.06047598e-04 3.80434787e-04 -9.99999462e-01]] Is it a good fit?2 x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff [[-2.77778132e-07 4.75306325e-10 8.22098410e-04] [ 4.75306325e-10 -2.77778132e-07 3.89069642e-04] [ 8.22098410e-04 3.89069642e-04 -9.99999462e-01]] Is it a good fit?1 x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff [[-3.25928641e-07 4.75306325e-10 9.65053368e-04] [ 4.75306325e-10 -3.25928641e-07 3.89726595e-04] [ 9.65053368e-04 3.89726595e-04 -9.99999462e-01]] Is it a good fit?1 x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff [[-3.89253374e-07 4.75306325e-10 1.15176103e-03] [ 4.75306325e-10 -3.89253374e-07 3.81414201e-04] [ 1.15176103e-03 3.81414201e-04 -9.99999462e-01]] Is it a good fit?2 x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff [[-4.61463753e-07 4.75306325e-10 1.36490746e-03] [ 4.75306325e-10 -4.61463753e-07 3.52580456e-04] [ 1.36490746e-03 3.52580456e-04 -9.99999462e-01]] Is it a good fit?2 x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff [[-5.47788559e-07 4.75306325e-10 1.61837109e-03] [ 4.75306325e-10 -5.47788559e-07 2.99621247e-04] [ 1.61837109e-03 2.99621247e-04 -9.99999462e-01]] Is it a good fit?3 x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff [[-1.64699627e-07 4.75306325e-10 5.20141144e-04] [ 4.75306325e-10 -1.64699627e-07 2.94129674e-04] [ 5.20141144e-04 2.94129674e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff [[-1.84993418e-07 4.75306325e-10 5.83879266e-04] [ 4.75306325e-10 -1.84993418e-07 2.95225429e-04] [ 5.83879266e-04 2.95225429e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff [[-2.07921155e-07 4.75306325e-10 6.56266515e-04] [ 4.75306325e-10 -2.07921155e-07 2.90163246e-04] [ 6.56266515e-04 2.90163246e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff [[-2.34890378e-07 4.75306325e-10 7.42205549e-04] [ 4.75306325e-10 -2.34890378e-07 2.78283250e-04] [ 7.42205549e-04 2.78283250e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff [[-2.64584114e-07 4.75306325e-10 8.37927431e-04] [ 4.75306325e-10 -2.64584114e-07 2.55953143e-04] [ 8.37927431e-04 2.55953143e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff [[-3.03670268e-07 4.75306325e-10 9.61501514e-04] [ 4.75306325e-10 -3.03670268e-07 2.22943889e-04] [ 9.61501514e-04 2.22943889e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff [[-3.31767296e-07 4.75306325e-10 1.05090498e-03] [ 4.75306325e-10 -3.31767296e-07 1.73365250e-04] [ 1.05090498e-03 1.73365250e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff [[-1.45522900e-07 4.75306325e-10 4.26517799e-04] [ 4.75306325e-10 -1.45522900e-07 2.63232521e-04] [ 4.26517799e-04 2.63232521e-04 -9.99999462e-01]] Is it a good fit?5 x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff [[-1.51287709e-07 4.75306325e-10 5.05660579e-04] [ 4.75306325e-10 -1.51287709e-07 2.39860394e-04] [ 5.05660579e-04 2.39860394e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff [[-1.65183726e-07 4.75306325e-10 5.53983800e-04] [ 4.75306325e-10 -1.65183726e-07 2.27513767e-04] [ 5.53983800e-04 2.27513767e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff [[-1.81567497e-07 4.75306325e-10 6.10196583e-04] [ 4.75306325e-10 -1.81567497e-07 2.11818232e-04] [ 6.10196583e-04 2.11818232e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff [[-1.97176154e-07 4.75306325e-10 6.65575727e-04] [ 4.75306325e-10 -1.97176154e-07 1.86834572e-04] [ 6.65575727e-04 1.86834572e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff [[-2.14258821e-07 4.75306325e-10 7.22344515e-04] [ 4.75306325e-10 -2.14258821e-07 1.57804687e-04] [ 7.22344515e-04 1.57804687e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff [[-1.65240971e-07 4.75306325e-10 5.30559597e-04] [ 4.75306325e-10 -1.65240971e-07 1.66345162e-04] [ 5.30559597e-04 1.66345162e-04 -9.99999462e-01]] Is it a good fit?5
plt.figure(figsize=(15, 10))
center_leds = center_cords_circle[3, 3]
y_unit_len = center_cords_circle[3, 2][1] - center_cords_circle[3, 3][1]
x_unit_len = center_cords_circle[4, 3][0] - center_cords_circle[3, 3][0]
center_leds, y_unit_len, x_unit_len
x_axis = [center_leds[0] - i*x_unit_len for i in range(-4, 4)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-4, 4)]
xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')
x_cords = [x if x != 0.0 else np.nan for x in center_cords_circle.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_circle.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range_circle.reshape(-1)]
plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()
count = 0
for x in range(12, 19):
for y in range(12, 19):
plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
count += 1
plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
axis_lens_pro = axis_lens_circle.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')
xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])
cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')
cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')
plt.show()
axis_lens_pro = axis_lens_circle.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
